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Algorithm is constructed which models single-file motion of particles 
interacting with each other and with the surroundings. As an example, 
we present the results of Brownian Dynamics simulations of the motion of 
cations moving through a short very narrow channel containing a device 
called "gate" , which may open and close the channel. 

1. Introduction 

Nanochannel transport and physical mechanisms of its regulation are 
among leading open problems in nanoscience. Its importance results from 
the fact that controlled and selective flow of matter through proteins in the 
cell membrane - achieved by active and passive channels - is one of most 
important biophysical processes in living cells. On the other hand, similar 
functions may be performed by synthetic nanopores which also can rectify 
the ionic currents [21 El H] and pump the ions against their concentration 
gradients |Sj and therefore may be used as simple models of biological (pro- 
tein) channels, and, on the other hand, may serve as devices for manipulat- 
ing the transport in the nanoscale. Therefore it is important to understand 
the conditions and properties of material transport inside the nanopore. 

The well-known and rather obvious property of the transport of material 
through very narrow pores is that the particles (ions, molecules ...) can pass 
through such channels in the form of single file only [H] . 

In the absence of noise (i. e., in standard Molecular Dynamics simula- 
tions) time increments 5t can be made arbitrarily small. This feature makes 
easy (in principle, at least) to keep all particles in prescribed unchanged or- 
der. In the Brownian Dynamics (BD) the action of random forces may result 
in arbitrarily high velocities and arbitrarily long jumps, time increment be- 
ing irrelevant in this respect. Therefore it is impossible to keep particles 
in still the same ordering by reducing time increments, the more that in 
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the presence of noise the time increments cannot be arbitrary IB] • Some 
additional procedures are needed. 

We present here the developed by us algorithm which models single-file 
motion of particles interacting with each other and with the surroundings, 
moving in a short very narrow channel containing a device called "gate", 
which may open and close the channel. To be specific, we shall discuss in this 
paper the electrostatic and hard-sphere interactions, though the formulas 
and the algorithms themselves can be easily adapted to any (sensible) form 
of interactions. 

2. The model 

We use the simplified model which does not take into account the details 
of the channel's structure. Full MD simulations of a K + -channel, including 
its molecular structure, water inside, all ions in the immediate vicinity, 
etc., requires use of total number of atoms in the simulation system above 
4 x 10 4 , and time-steps 0.2 fs (HI Ell- Such simulations have also some other 
drawbacks 

Little is known about the details of the gating mechanism, the more 
that the motions of dangling ends ^2] i n synthetic pores are probably quite 
different from the motions of the subunits of proteins constituting the bi- 
ological channels. Therefore, without entering into details of equations of 
motion for the channel's walls, we model the gating process by introducing 
inside the channel the artificial device called "gate" which can either allow 
or prevent the flow of particles through the channel. 

The main assumptions are: 

(i) We simulate the motions of the particles inside the simulation zone 
(SZ) of the lenght L, narrow enough to force the particles inside SZ to move 
in the single-file order. Knowledge of the detailed shape (e. g. cylinder, 
cone, hour-glass) is not necessary from this point of view. Regions outside 
are treated as reservoirs for particles both outcoming from and ingoing into 
SZ. 

(ii) We neglect the motions in radial directions, and describe the parti- 
cles as moving along the z-axis of the SZ only (quasi-onedimensional mo- 
tion). However, the physical system (electrostatic interactions, etc.) re- 
mains three-dimensional. 

(iii) The opening and closing of the channel (so-called gating process) 
is modeled by the presence of the charged "gate" located inside SZ. The 
state of the gate is determined by its Brownian motion (Wiener process of 
intensity Qb), and by electrostatic interactions with the ions inside SZ and 
with external electric field. The gate opens when the net force exceeds some 
threshold value, and closes otherwise. Minimal approach distance between 
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particle and gate is d cg . 

(iv) The real channels exhibiting the flicker noise are asymmetric and 
charged. We model these properties by the mentioned above gate, and by 
additional charges located outside SZ. 

(v) Water molecules are not modeled explicitly but are described elec- 
trostatically by an effective dielectric constant and as the source of friction 
and noise - as is frequently done [P5] . 

No periodic boundary conditions are imposed. Instead, in our simula- 
tions we assumed (when other rules are satisfied) that 

(i) Particles can leave and enter the simulation zone (SZ) through both 
apertures. 

(ii) Particle leaves the simulation zone (and can be counted to the current 
balance at the given aperture) when its center-of-mass position is smaller 
than the lower threshold, or greater then the higher threshold. In our case 
we accepted as thresholds the particle diameter d c and SZ length L minus 
d c . 

(iii) Single-file assumption implies that when one particle leaves the sim- 
ulation zone, another cannot enter through the same aperture in the same 
time (i. e., during the same time-step). 

(iv) When rule (iii) allows, particle may enter SZ when nearest particle 
is farther that the prescribed smallest distance. In our case the smallest 
distance is d c + e (e = 0.00001 nm). 

(v) Particles enter SZ with prescribed finite probabilities -P(O) and P(L), 
which may be different for different apertures (i. e. at x = and x = L). 
The probabilities of entrance simulate concentrations outside SZ - the lower 
concentration, the lower probability. 

In our simulations we assumed that (when other rules are satisfied) 
during one time-step only one particle may enter the SZ through a given 
entrance, and, when entering, that it is located at the distance d c from the 
aperture. This rule can be changed. 

The Langevin-type equations of motion for the particles (cations) mov- 
ing along the channel reads: 

mvi = -jm + Ri(Zi) + Fi(zi) , 

Zi = Vi, (1) 

where Vi is the velocity of i-th. ion, Zi - the position, m; - the mass, ji - 
the friction coefficient, Fi(zi) - sum of deterministic forces, and R%{zi) — the 
random force. 

The gate is charged to prescribed value q g = Z g e, where Z g is the valence 
and can be in two states: open and closed, respectively. In our simulation 
important is the absolute value of the force F g acting on the gate. We 
assume F g to be sum of deterministic and random forces described below. 
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The deterministic force Fi(zt) experienced by the cations and the gate 
consist of the applied external force (voltage), and the internal Coulomb 
force from other charges. The Coulomb interaction between two ions is 
modified by the addition of a short-range repulsive 1/r 10 force, where r is 
the ion-ion distance [S]. 

The random force Ri acting on ions is assumed to be the thermal noise 
represented by the Gaussian white noise. On the other hand the random 
force experienced by the gate R g is given by the Wiener process (gate's 
Brownian motion) R g = J2i Ri- 

In the Brownian Dynamics calculations, St should be of the order of 771/7 
|15l [7| • Using the Euler scheme 

mv(t)+^v(t)=F(t) -> m ^- + 5t ) - ^ + jv(t) = F(t) (2) 

ot 

would lead to obviously wrong result: v (t + St) = F(t). Therefore we use 
the following scheme of discretization: 



v(t + St)—v(t) 7. , ... 
j t ±l + ±[v(t + St)+v(t)} = F(t) 

z(t + St) - z(t)) 

= V (t + St). (3) 

This computational scheme is similar, though not identical, with that de- 
scribed recently in ref. [TIE]- The "forward evaluation" (Eq. 3) has stability 
and accuracy implications, and ^1] suggest using it for each extrapolative 
force calculations. 



3. Numerical results 

The length of the simulation zone is L = 10 nm. This corresponds to 
the real length of biological channels, and - roughly - to the length of the 
narrow part of the synthetic channel reported in 

The net flow of particles through the channel (simulation zone) was cal- 
culated either by keeping the balance of particles entering and leaving both 
apertures, or by counting the particles passing the gate in both directions. 
Both procedures lead to the same results. 

Initial values of velocities of particles were drawn from the Maxwell dis- 
tribution with the variance ksT '/m c . The results are insensitive on the exact 
values of temperature and mass within rather wide range of temperatures 
and masses. 

A list of the parameters used in the BD simulations is given below: 
Temperature: T = 298 K and k B T = 4, 12 x 1CT 21 J, 
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Mass: m c = 6.5 x 10~ 26 kg, Friction constant: "f c = 2.08 x 1CT 12 kg/s, 
Dielectric constant: e w = 81, Voltage: U = 1.77 x 1CT 2 V 
Ion diameter: d c = 0.266 x 10 -9 m, Ion-gate min. distance: d cg = 2.5c£ c , 
Valences: Z c = +1, Z g = —50, 
Intensity of short-range force: F% R = 444 x 10" 9 N, 
Intensity of noise: Q c = 0.47 x 10~ 9 N, Q g = OMQi. 
Intensity of the gate's noise Q g is different from cations' one Q c (and 
is taken as a free parameter) due to the difference of masses, and also due 
to a kind of "stiffness" of (or hindrances in) the motions of channel's walls 
constituents. In all simulations first 10 6 steps were rejected. The power 
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Fig. 1. Power spectra S(f) of the stochastic series of subsequent values of the 
net number of cations m n leaving the simulation zone. Red: S(f) for 7 different 
realizations of the intrinsic noises, the same values of all parameters in every series: 
e = 81, m c = 6.5 x 1(T 26 kg, U = 17.7 * 10~ 3 V, St = 31 x lCT 15 s, Q g = 0.01Q C , 
F° R = 444 x 1CT 9 N, gate thresholds = ±1100 x 10~ 12 N. Blue: S(f) with the same 
realizations of the intrinsic noises, for 7 different values of all parameters in every 
series. In every series only one parameter is changed: e = 0.93e°, m = 0.77m , 
U = 1.33J7 , St = O.SSt , Q g = 0.75Q°, F S r = 0.6F% R , where p° denotes the value 
of the given parameter from the panel A. 



spectrum was calculated from runs of length W 7 5t. The power spectrum of 
the series {mi, . . . mjv} is 



S <" = N 



N 

£ m n e- 2m f n 

n=l 



(4) 



where m n denotes either the net number of particles leaving SZ during the 
n-th step (then m n can be either positive, zero, or negative), the number 
of particles inside SZ at the end of the n-th step (m n = N p > 0), or 
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the state of the gate during the n-th step (then m n = {0,1}). All these 
power spectra are dimensionless. There are data that suggest that inside 
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Fig. 2. Power spectra S(f) of the stochastic series of subsequent values of the 
number of cations N p inside the simulation zone. Notation and values of parameters 
the same as in Fig. 1. 

very narrow pores the physical properties of aqueous solutions, such as 
dielectric constant, density, diffusion coefficient, viscosity, solvatation of ions 
(i. e., their effective diameters), etc. may differ from their bulk values |17j . 
Therefore we checked how the changes of such parameters influence our 
model. We found that the quantitative changes of calculated values of net 
currents and of frequency spectra resulting from reasonable variations of 
these parameters are within the limits of quantitative differences resulting 
from different realizations of the noise. The results are shown in Figs. 1-3. 
These observations suggest robustness of the model. On the other hand, 
the model is sensitive with respect to the changes of relative strength of 
random and deterministic forces - decrease of the dielectric constant with 
noise unchanged, or increase of noise with electrostatic forces unchanged 
changed significantly the results. E. g., either too strong gate noise or too 
strong electrostatic force (i. e., low dielectric constant) dampen the flicker 
noise. 

When the single-file limitations are removed, all the power spectra shown 
in Figs. 1-3 become S(f) ~ i -e, the corresponding processes behave 

like the Wiener process. 

4. Appendix 

Here we present the codes for the single-file motion. The codes for 
entrances and exits of particles, for the number of particles located to the 
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Fig. 3. Power spectra S(f) of the stochastic series of subsequent values of the state 
of the gate. Notation and values of parameters the same as in Fig. 1. 

left of the gate, as well as the codes for the determination of the state of 
the gate (open or closed), and for the equations of motion are standard and 
will not be reproduced here. 

Single-file procedures are based on the fact that the given particle (cation) 
i cannot move farther that its neighbours i — 1 and i + 1, which in turn are 
limited by their neighbours, i and i — 2 or i + 2, etc. Therefore their po- 
sitions need to be recalculated. In the simplest version, it is assumed that 
particles meet at the middle of their former positions. In the better versions 
such a pair of particles meets at the position calculated from their former 
positions and from their new velocities. On the other hand, the particles 
retain their velocities until a given pair meets, then they collide and - in 
the simplest version - exchange their velocities (behave as hard spheres). 
Again, it is possible to refine this simplest procedure. Because the results 
of the above-described procedure depend on whether the recalculations are 
done "up" or "down", i.e., from particle number 1 to N, or from N to 1, 
both reorderings are realized independently, their results are averaged, and 
the whole scheme is iterated until self-consistency is attained. 

Before using the SFM-codes below, one needs to supply the values of 
entries of three main arrays: ZK[Nkmax], VK[Nkmax], in which the po- 
sitions and velocities of particles inside the simulation zone are stored, 
ZK0[Nkmax] in which former positions are remembered, and two auxiliary 
ones: ZKG[Nkmax] and ZKH[Nkmax] for storing intermediate data. It is 
needless to say that these arrays should be declared as external variables. 

Nkmax denotes here the maximal, Nk (in the codes) - the actual number 
of particles inside the simulation zone. 

We present here separate single- file codes for open and for closed channel. 
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Before calling the S-F code for a closed channel, one needs to calculate the 
number of particles located to the left of the gate, denoted in the codes as 
Nkgl. 

// single-file ordering: open channel 
nrep = ; repeat = 1 ; 
while (repeat == 1 && Nk > 1) 

{ 

repeat = 0; nrep++; 

orderlow(l,Nk,dcc) ; //pairs (1,2) , . . . (Nk-l,Nk) 
orderup(0,Nk-l,Nk,dcc) ; //pairs (Nk,Nk-l) , . . . (2, 1) 
ave(l,Nk) ; 
subst(Nk) ; 

} 

// single-file ordering: closed channel 
while (repeat ==1) 

{ 

repeat = 0; nrep++; 

// ordering to the right of the gate: 
if (Nkgl < Nk) 

{ 

ic = ordergr(bp,Nkgl,Nkgl+l,Nk,dcc) ; // at the gate 

if(ic < Nk) // remaining particles 

{ 

orderlow(ic,Nk,dcc) ; // pairs ic, ic+1) , . . . (Nk-1 ,Nk) 

orderup(0,Nk-ic,Nk,dcc) ; //pairs Nk, Nk-1) (ic+1 , 

ave (ic ,Nk) ; 

} 

// ordering to the left of the gate: 

} 

ic = ordergl (bl, Nkgl, 0, Nkgl, dec) ; // at the gate 

if(ic > 1) // remaining particles 

{ 

orderlow(l,ic,dcc) ; // pairs (1 ,2) , . . . (ic-1 , ic) 
orderup(0,ic-l,ic,dcc) ; // pairs (ic-1 , ic) , . . . , (2 , 1) 

ave (1 , ic) ; 

} 

subst(Nk) ; 

} 

void orderlow(int m, int N, double d) 

{ 
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int i ; 

f or(i=m; i<N; i++) 

{ 

if(ZK[i] > ZK[i+l] - d) 

{ 

ZKH [i] = 0.5*(ZK0[i] + ZKO[i+l] - d) ; 
ZKH[i+l] = ZKH [i] + d + 0.00001; 
VKH [i] = VK[i+l]; VKH[i+l] = VK[i] ; 
repeat = 1 ; 

} 

} 

return; 

} 

void orderup(int m, int N, int M, double d) 

{ 

int i , j ; 
for(j=m; j<N; 
{ i = M - j; 

if(ZK[i] < ZK[i-l] + d) 

{ 

ZKG [i] = 0.5*(ZK0[i] + ZK0[i-l] + d) ; 
ZKG[i-l] = ZKG [i] - d - 0.00001; 
VKG [i] = VK[i-l]; VKG[i-l] = VK[i] ; 
repeat = 1 ; 

} 

} 

return ; 

} 

int ordergl (double b, int Nkgl, int m, int N, double 

{ 

int i, j , ii; 
double a, c; 

be = b; ii = Nkgl; // be = gate + ded 
for(j=m; j<N; 

{ 

i = Nkgl - j ; 
if(ZK[i] > be) 

{ 

ZK[i] = be; ZKG [i] = be; ZKH [i] = be; 
be -= d; ii = i; 
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if(i == Nkgl) VK[i] = -VK[i]; else 

{ 

a = VK[i] ; c = VK[i+l] ; 
if(fabs(a) > fabs(c)) 

{ 

VK[i] = a + c; VK[i+l] = 0; 

} 

else 

{ 

VK[i+l] = a + c; VK[i] = 0; 

} 

} 

repeat = 1; 

} 

else break; 

} 

return ii; 

} 

int ordergr (double b, int Nkgl, int m, int N, double 

{ 

int i, ii; 

double a, c; 

be = b; ii = Nkgl + 1; 

f or(i=m; i<=N; i++) 

{ 

if(ZK[i] < be) 

{ 

ZK[i] = be; ZKG [i] = be; ZKH [i] = be; 
be += d; ii = i; 

if(i == Nkgl+1) VK[i] = -VK[i]; else 

{ 

a = VK[i] ; c = VK[i-l] ; 
if(fabs(a) > fabs(c)) 

{ 

VK[i] = a + c; VK[i-l] = 0; 

} 

else 

{ 

VK[i-l] = a + c; VK[i] = 0; 

} 

} 
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} 

else break; 

} // i 
return ii; 

} 

void ave(int m, int N) 
{ 

int i ; 

f or(i=m; i<=N; i++) 
{ 

ZK[i] = (ZKG [i] + ZKH[i])/2.0; 
VK[i] = (VKG [i] + VKH[i])/2.0; 

} 

return; 

} 

void subst(int N) 

{ 

int i ; 

for(i=l;i<=N;i++) 

{ 

ZKG [i] = ZK[i] ; ZKH [i] 
VKG [i] = VK[i] ; VKH [i] 

} 

return; 

} 

For simplicity, we give here, in the functions orderlow and orderup 
(lines ZKH [i] = 0.5(ZK0[i] + ZK0[i+l] - d) , ZKG [i] = 0.5(ZK0[i] + 
ZK0[i-l] + d)) the simplest form of the recalculations of the correct po- 
sitions of pairs of particles as contact positions of the pair in the middle 
of their former positions. These positions can be determined with better 
accuracy by taking into account particles' velocities, and their equations of 
motion as well. The appropriate codes are obvious and are not reproduced 
here. 

These codes can be written in a more compact way. The form presented 
here is - in our opinion - better legible and more self-explanatory. 



= ZK[i] ; 
= VK[i] ; 
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